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Scattering  of  elastic  waves  by  a  spherical  inclusion  -  1.  Theory  and 

numerical  results 

Valeri  A.  Korneev  and  Lane  R.  Johnson 

Department  of  Geology  and  Geophysics,  University  of  California, 
and  Center  for  Computational  Seismology,  Lawrence  Berkeley  Laboratory, 
Berkeley,  California  94720 

A  complete  and  exact  solution  for  the  problem  of  an  incident  P  wave  scattered 
by  an  elastic  spherical  inclusion  is  presented  and  described.  The  solution  can  be 
obtained  from  either  analytical  formulas  or  stable  numerical  procedures.  A  method  of 
estimating  the  number  of  terms  that  must  be  retained  in  the  harmonic  series  in  order 
to  achieve  a  specified  accuracy  is  given.  Lite  results  are  investigated  by  calculating 
synthetic  seismograms,  scattering  diagnuns,  ;uid  scattering  cross  sections  for  a  broad 
frequency  band  and  for  both  low-velocity  ;ind  high-velocity  inclusions.  The  fields 
within  the  shadow  zone  are  formed  primarily  from  three  different  types  of  waves,  P 
waves  transmitted  through  the  sphere,  P  waves  diffracted  around  the  sphere,  and  S 
waves  converted  at  the  boundary  of  the  sphere.  The  relative  contribution  from  these 
different  waves  depends  upon  the  distance  of  the  observation  point  from  the  sphere. 
Key  words:  elastic  waves,  sphere,  scattering,  diffraction 


1.  Introduction 

The  fact  that  the  earth  is  not  a  homogeneous  body  has  forced  seismologists  to  consider  scattered 
waves  in  their  attempts  to  explain  some  of  the  features  they  observe  on  seismograms.  Elastic  wave 
scattering  has  been  called  upon  to  explain  a  variety  of  different  phenomena  which  are  routinely 
observed.  These  include  the  phase  and  amplitude  fluctuations  of  waves  arriving  at  a  seismic  array  (Aki. 
1973),  the  precursors  to  PKIKP  (Haddon  and  Cleary,  1974),  and  the  codas  of  local  earthquakes  (Aki, 
1969).  in  addition,  the  attenuation  of  seismic  waves  is  usually  interpreted  as  being  a  combination  of 
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intrinsic  absorption  and  scattering  (Aki,  1980). 


The  general  problem  of  scattering  of  elastic  waves  by  a  heterogeneity  within  the  earth  is  a 
difficult  one  and  analytical  solutions  are  known  for  only  a  few  special  cases,  and  even  in  these  instances 
u>e  solutions  are  complicated  and  laborious  to  calculate.  Thus,  most  attempts  at  interpreting  scattered 
seismic  waves  have  relied  upon  approximate  treatments  of  the  scattering  theory.  A  number  of  approxi¬ 
mations  are  possible  and  are  related  to  such  parameters  as  the  size  of  the  scatterer,  the  shape  of  the 
scatterer,  the  distance  of  the  observation  point  from  the  scatterer,  the  magnitude  of  the  heterogeneity, 
and  the  number  of  times  the  wave  has  been  scattered.  Assessing  the  validity  of  these  approximations  is 
not  a  simple  matter  and  has  been  addressed  for  only  some  of  the  approximations  (for  example,  Hudson 
and  Heritage,  1981).  The  matter  is  further  complicated  by  the  fact  that  more  than  one  approximation 
may  be  involved  in  the  stune  problem,  and  it  is  not  always  obvious  that  the  different  approximations  are 
coasisient. 

The  best  method  of  checking  the  validity  of  the  approximations  which  are  made  in  scattering 
problems  is  to  compare  them  with  exact  analytical  solutions.  In  this  paper  we  develop  and  discuss  the 
properties  of  one  such  solution,  the  scattering  of  a  plane  P  wave  by  a  spherical  inclusion.  The  method 
of  obtaining  the  analytical  solution  is  outlined,  the  numerical  methods  used  in  calculating  the  solution 
are  described,  and  some  of  the  important  features  of  the  solution  are  described  and  discussed.  The 
results  contained  in  this  paper  ;ire  the  suit  Ling  point  for  a  companion  paper  (Korneev  and  Johnson, 
1992),  which  makes  detailed  comptirisons  between  various  approximations  and  the  exact  solution. 

In  seUing  up  the  theoretical  problem  for  elastic  wave  scattering,  one  must  choose  a  model  which 
is  complicated  enough  so  that  it  resembles  situations  found  in  the  earth  and  at  the  same  time  simple 
enough  to  allow  tractable  solutions  of  the  mathematical  equations.  We  have  chosen  to  model  a  local 
heterogeneity  as  a  homogeneous  elastic  sphere  surrounded  by  a  homogeneous  elastic  medium.  The 
elastic  constants  and  density  of  the  spherical  inclusion  can  be  ;irbitr;irily  different  from  those  of  the  sur¬ 
rounding  medium  The  sphere  can  also  be  Idled  with  a  lluid.  The  sphere  is  one  of  the  few  objects  for 
which  the  scattering  problem  can  be  solved  exactly.  It  also  has  the  desirable  property  of  being  describ- 
able  by  a  minimum  number  or  parameters,  which  makes  the  interpretation  of  the  analytical  and 
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numerical  results  relatively  simple.  Furthermore,  scattering  by  a  sphere  represents  a  canonical  problem 
for  a  more  extended  class  of  objects  with  relatively  simple  and  smooth  boundaries,  and  thus  the  results 
presented  in  this  paper  contribute  to  the  general  understanding  of  scattering  from  this  class  of  hetero¬ 
geneities. 

2.  Theory 

The  analytical  treatment  of  scattering  of  elastic  waves  by  a  spherical  inclusion  has  a  long  and  rich 
history,  dating  back  to  the  classical  papers  by  Clcbsch  (1863)  and  Rayleigh  (1871).  Much  of  the  earlier 
work  was  restricted  to  scalar  (acoustic)  waves,  hollow  inclusions,  or  rigid  inclusions.  Much  of  the  more 
recent  work  on  the  scatterering  of  elastic  waves  by  a  spherical  inclusion  has  been  associated  with  the 
names  of  G.  I.  Petrashen  in  the  USSR  and  R.  Truell  in  the  USA.  As  early  as  1944  Petrashen  was 
studying  the  scattering  problem  for  an  incident  P  wave,  where  he  made  use  of  a  natural  spherical  vector 
system  which  he  h;td  developed  earlier  for  the  purposes  of  quantum  mechanical  problems.  The  later 
studies  by  Petrashen  and  his  students  were  concerned  mostly  with  obtaining  asymptotic  representations 
for  different  kinds  of  regular  and  diffracted  waves  formed  by  a  sphere  (Petrashen,  1946,  1950a  1950b, 
1953;  Buldurev  and  Molotkov,  1958),  such  as  head  interference  waves  for  the  scalar  case  (Buldyrev  and 
Molotkov,  1960),  failure  waves  for  the  scalar  case,  surface  waves  on  an  isolated  elastic  sphere  (Gelchin- 
skij,  1958),  and  complete  solutions  for  an  elastic  sphere  (Korneev,  1983;  Korneev  and  Petrashen,  1987). 
Related  work  includes  the  theoretical  and  experimental  studies  of  acoustic  waves  incident  upon  an  elas¬ 
tic  sphere  by  Nigul  et  al.  (1974).  Independently,  Truell  and  his  coworkers  (Ying  and  Truell,  1956;  Ein- 
spruch,  Witterholt,  and  Truell,  1960;  Truell,  Elbaum,  and  Chick,  1969)  also  made  important  contribu¬ 
tions  to  the  solution  of  this  problem.  The  basic  equations  that  are  obtained  when  displacement  poten¬ 
tials  are  used  along  with  summaries  of  some  of  the  work  mentioned  above  can  be  found  in  such  books 
as  Morse  and  Feshbach  (1953)  aid  Pao  and  Mow  (1971). 

The  basic  method  followed  in  the  above  papers  is  to  write  the  solutions  inside  and  outside  the 
inclusion  in  terms  of  appropriate  eigenfunctions  of  the  differential  equations  and  then  couple  these  solu¬ 
tions  by  matching  the  boundary  conditions  on  the  surface  of  the  inclusion.  It  is  also  possible  to 
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formulate  this  problem  as  an  integral  equation  by  making  use  of  the  elastodynamic  Green  function  for  a 
homogeneous  medium.  Further  information  on  this  approach  can  be  found  in  Miles  (1960),  Haddon 
and  Cleary  (1974),  Varatharajulu  and  Pao  (1976),  Wateiman  (1976),  Hudson  (1977),  and  Gubematis  et 
al.  (1977a,  1977b). 

In  this  section  we  oudine  a  method  of  obtaining  the  general  solution  for  the  scattering  of  a  plane 
P  wave  by  a  spherical  inclusion.  Although  a  solution  to  this  problem  can  be  found  in  several  of  the 
papers  listed  in  the  previous  paragraphs,  we  repeat  it  here  for  completeness  and  also  to  introduce  a  set 
of  spherical  basis  'ectors  which  are  used  in  our  solution  but  are  not  common  in  the  English  literature. 
These  basis  vectors  were  developed  by  Petrashen  and  have  properties  which  are  convenient  for  solving 
vector  function  problems  with  spherical  symmetry. 

Consider  the  configuration  shown  in  Figure  1.  Joint  Cartesian  {jc  ,y,z)  and  spherical  (r,0,<|>) 
coordinate  systems  will  be  used.  Centered  about  the  common  origin  of  these  coordinate  systems  is  a 
sphere  of  radius  r~R  The  volume  within  this  sphere  will  be  denoted  by  the  index  v  =  1 ,  while  the 
volume  outside  will  be  denoted  by  v  =  2.  The  materials  inside  and  outside  the  sphere  will  in  general 
be  different,  and  the  properties  of  these  materials  are  completely  described  by  the  Ltune  parameters  and 
density 

K  =  const.  .  pv  =  const.  ,  pv  =  const.  ,  (v  =  1,2)  (2.1) 

Incident  from  medium  v  -  2  is  a  harmonic  disturbance  with  a  displacement  field  given  by 

fJ„  =  U„  U.v (2.2) 

The  interaction  of  this  incident  wave  with  the  sphere  gives  rise  to  additional  displacement  tields  both 
inside  and  outside  the  sphere,  and  these  are  denoted  by 

Uv  =  Uv  (x,y.z)e,ov  .  (v  =  1,2)  (2.3) 

These  ttddiliontd  disturbtuiccs  will  be  referred  to  as  the  scattered  fields.  Associated  with  each  field  is  a 
stress  tensor.  We  will  only  need  the  traction  on  spherical  surfaces  and  this  is  given  by 

dUv 

t,!v)(Uv)  =  Xv  V-  Uv  f  +  2pv  — —  +  p^r  X  Vx  Uv)  (2.4) 

or 
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where  F  =  rtr  is  the  unit  radius-vector.  We  also  denote  the  velocities  and  slownesses  of  both  compres- 
sional  ;uul  shear  waves  by 


VXv+2fiv 
Pv 


=  Otv'1 


V:  ■  *■ 


The  incident  field  of  equation  (2.2)  and  the  scattered  fields  of  equation  (2.3)  must  all  satisfy  the 
elastodynamic  equations  of  motion  in  mediums  v  =  1,2  , 

(Xv+2ftv)  V’Uv  -  |iv  VxVxl!v  +  pvurUv  =  0  (2.6) 

Taken  together  the  fields  must  satisfy  the  boundary  conditions  on  the  surface  of  the  sphere,  which  tire 
that  the  displacement  ;ind  traction  should  be  continuous.  Thus  we  require  that 


U,  =  U„  +  U:  and  t,(l)  U,  =  t<2!ju(1  +  U2 


where  all  of  these  equations  are  evaluated  at  r  =  R .  We  also  require  that  the  scattered  fields  remain 
finite  within  the  sphere  and  satisfy  a  radiation  condition  at  targe  distances  from  the  sphere 

U2  =  r-*'  (  r  -»  oo)  (2.8) 

where  k  =  (i)/V^2)  or  k  =  co/V,l2).  This  is  a  well-posed  problem  in  that,  given  the  incident  wave,  the 
boundary  conditions  are  sufficient  to  solve  for  the  scattered  fields  and  thus  arrive  at  a  unique  solution  to 
the  problem. 

Now  let  us  consider  the  special  case  of  an  incident  wave  which  is  a  plane  harmonic  P  (compres- 
sional)  wave  propagating  in  the  direction  of  the  positive  z  axis.  In  this  case  we  have 

..  -Mini-,:  „ 

u<>  =  C  -  Z  (2.9) 

where  z  is  a  unit  vector  in  the  z  direction  We  construct  the  solutions  by  introducing  the  system  of 
spherical  vectors  developed  by  Petrashen  (1945.  1949) 


y;„  -Y£,fe.o)  =  rx  VK,„,(e,0) 


Yu  -  Y£<0.0)  =  (Ml)  r  K,m(0.< ft)  -  r  VKu(0.O> 


Yu  =  Yu  (0,0)  =  /  f  Y im  (0.0)  +  r  VK,„,(0.O) 
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with  the  usual  definitions  for  the  spherical  harmonic  functions 


=  em*  PP(cosQ)  .  t  >  0  ,  (-/  <  m  <  l) 


This  system,  whose  main  features  are  outlined  in  Appendix  A,  leads  to  particularly  simple  equations  in 
this  problem  of  spherical  symmetry,  Using  the  system  of  equation  (2.10),  an  arbitrary  vector  function 
U  can  be  represented  in  the  form 


U  =  2  j  <(r)  M  +  V£(r)  ML  +yL(r)  Yfm 


where  the  v |f/,„(r)  ;ue  unspecified  radial  functions  at  this  point 


The  particular  incident  wave  of  equation  (2.9)  has  the  representation 


U<>  =  2j  y/+i(toa2r)  Y,{, -;,.|((oa2r)  Y,o  •< 
/20 


where  the  jk(x)  are  spherical  Bessel  functions.  The  displacement  fields  of  the  scattered  waves  U,  and 
ll2  Irom  equation  (2..))  could  also  l>c  rcpicseutcd  in  terms  ol  a  series  of  the  sphciical  vcitois  ol  equa¬ 
tion  (2.10).  Making  use  of  the  structure  of  the  series  in  equation  (2.12)  for  the  incident  wave  and  the 
orthogonal  properties  of  the  spherical  vectors  (see  Appendix  A),  it  is  possible  to  show  that  U,  and  U2 
can  be  represented  by 


U|  =  (t/'V^ifoxxi/-)  1-  //V'V/+ j04V)  Y/o 


+  — ci/"7/_i(coot,r )  +  (/  + 


-I  —  (/♦!) 

l)/r/(lV,-l((uPir)J  Y,',  -e  2 


If)  tiia>liitt((im,r)  •  //;/2,/i,h((o |Ur)  Y,*, 


r  ,1  -i— ir+i 

+  |-«/‘>/i,_|(<i)a,r)  +  (/  +  D/V2lh/.|(a)(5ir)j  Mf„  \e  2 

j 
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where  a/v)  and  b/v’  are  now  unknown  scalars  and  the  hk(x)  are  spherical  Hankel  functions  of  the 
second  kind.  Note  that  the  spherical  Bessel  functions  jt U )  of  equation  (2. 13)  are  finite  throughout  the 
inclusion  and  the  spherical  Hankel  functions  hk(x)  of  equation  (2.14)  have  asymptotic  representations 
for  large  argument  that  satisfy  the  condition  in  equation  (2.8).  Also  note  that  with  the  basis  vectors 
being  used  here,  the  displacements  in  equations  (2.13)  and  (2.14)  separate  into  compressional  and  shear 
fields,  so  that  either  equation  can  be  put  in  the  form 

Uv  =  U''  +  lly  (2.15) 

with  the  a/v)  associated  with  the  P  field  U£  and  the  fr/v)  associated  with  the  S  field  Uv,  and  with  the 
following  conditions  satisfied 

Vx  U f!  =  0  .  V  Uv*  =  0  .  (2.16) 


To  detenu  mate  the  coefficients  a/v,and  b-[v>  in  equations  (2.13)  and  (2.14)  we  have  to  satisfy  the 
boundary  conditions  of  equation  (2.7).  Making  the  necessary  substitutions,  evaluating  the  expressions 
for  r  =  R ,  and  using  the  orthogonal  properties  of  the  spherical  vectors,  one  arrives  at  a  separate  set  of 
linear  equations  for  each  value  of  / .  When  /  =  0  one  can  take  advantage  of  the  fact  that  Y,*,  =  0  to  get 
the  abbreviated  set  of  equations 


W 


0 


(2.17) 


where 


W„ 


i(4i)  -  M^>> 
ko,  -  d  . 

v  *  ~  J 


(2.18) 


7ti(^v)  -  4 j,(%v) 
Yv 


%  A„(^:)  -  4h  ,(^2) 
Y i 


^  =  ct vw/? 


T)v  =  Pv0)« 

P1P2 
*  "  P:Pi 


(2.19) 
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When  l  >  la  full  set  of  four  equations  with  four  unknowns  is  obtained  of  the  form 


W,  F,  =  Ft 


where  Ft  and  Ft  are  the  column  matrices 
F,  =  [a/11,  ft/1*,  a,m ,  ft/2^ 


and 


Ft  -  [ft+i(^2)  .  ~jt- 1(^2)  •  C*  ,  C2 


and  the  matrix  Wt  is  given  by 


(2,20) 


(2.21) 


(2.22) 


Ut*  M\) 

-ift;+i(T)2) 

-jl- l(4l) 

(/+l)7/.i(tii) 

hi- 1(^2) 

-(f+1)  ft/-|(ri2) 

kC,+ 

K/  D,+ 

-c2+ 

-/  d2+ 

KCf 

-K(/  +  l)Df 

-c2- 

(/+dd2- 

The  following  notations  tire  used  here 


c;l  yL  J(/+2)  ;,♦,(&,)  1 

c;J  =  ^7  y,(4v)- 2  l</-d ;,->($.)  J 


c; 

c; 


£>v+ 

0; 


o: 

0; 


y 


[  u+2)  ft/+1(ul 

—  ftt(^v)  -  2 
Yv 

[  J 

(/+2)ri+|0lv)] 

»lv  ririlv)-  21 

(/ — 1 )  y/_i(riv)  J 

[  1 

(1+2)  ft/  +  i(T)v)  1 

Pv  ft/(qv)  -  2  ' 

(/-l)  /l;-|(qv)  J 

(2.23) 


(t  can  be  si.,.wn  dial  the  detenninanfs  of  equations  (2.  IX)  and  (2.23)  arc  always  different  from 
zero,  so  in  principle  11  is  always  possible  to  solve  the  systems  of  equations  (2.17)  and  (2.20)  for  the 
unknown  coefficients  «/V|  and  ft/V|  for  any  /  However,  because  of  the  asymptotic  properties  of  the 


X 


mainx  elements,  numerical  difficulties  can  he  encountered  when  solving  these  systems  of  equations. 
This  will  be  discussed  in  the  next  section. 

Finally,  note  that  the  coefficients  n/V)  and  b/V|  for  any  /  are  completely  defined  by  the  member 
of  the  series  in  equation  (2.12)  for  the  incident  wave  with  the  same  /.  That  is,  there  is  no  coupling 
between  harmonics  associated  with  different  values  of  l .  This  means  that  a  scattering  problem  for  the 
sphere  with  an  incident  wave  represented  by  any  member  of  equation  (2.12)  with  index  l  could  be  con¬ 
sidered,  and  the  results  would  be  the  corresponding  members  of  equations  (2.17)  and  (2.20)  with  the 
same  index  / .  To  illustrate  this  point  we  will  show  how  the  solution  for  an  incident  wave  generated  by 
a  point  pressure  source  can  be  obtained  by  a  slight  modification  of  the  solution  for  a  plane  wave  source 
developed  above.  Consider  a  point  pressure  source  located  at  the  point  R„  =  (Z(>,(),0)  where  Z0  >  R . 
This  will  generate  a  displacement  field  ll0  of  a  spherical  P-wave 

r)  Yfn  -  r>  Y,oj  h,(kpZ„,  (2.24) 

Each  member  of  the  series  in  equation  (2.24)  differs  from  that  in  equation  (2.12)  by  the  coefficient 

Ci  =  -i.:‘  e  -  hi  (fpZo) 

Thus  the  solution  for  the  point  pressure  source  is  simply  obtained  by  multiplying  the  members  of  equa¬ 
tions  (2.13)  and  (2.14)  by  q .  The  equations  (2.17)  and  (2.20)  are  unchanged. 


=  -v- 


-dt  lr-Rnl 


Ir-Ro 


= 


Vi«(k. 


IX) 


3.  Numerical  considerations 

To  calculate  the  scattered  fields  LI,  and  U2  using  the  scries  in  equations  (2.13)  and  (2.14)  we 
have  to  solve  the  line;ir  systems  of  equations  (2.17)  and  (2.20)  for  all  /  >  0  .  These  systems  of  equa¬ 
tions  are  generally  solved  by  numerical  methcxls,  which  means  that  the  numerical  stability  of  the  calcu¬ 
lations  are  a  concern.  The  coefficients  of  the  matrices  W,  in  these  systems  include  spherical  Bessel 
and  Hankel  functions,  whose  features  for  /  >  1  are  characterized  by  the  formulas 

1  n  1  -'<*  -  t(/  + 1,1 

J,(Z)  ~  -  mv|  ;  -  -(/  +  1)  1  .  h,(z)  =  -  e  2  (3.1) 

z 
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if  z  »  l,  and 


ji(z)  = 


2z^y 


2l 


-<'4> 


hi(z) 


= 

zvy  [  < 


'4 


(3.2) 


if  z  «  /  ,  where 


y  = 


,  +  i 


Because  of  different  ratios  between  /  and  the  arguments  4i>  £2.  0i.  and  T|2  of  the  functions  j,  and  />/ 
that  appear  in  equation  (2.23),  the  orders  of  the  different  columns  of  this  matrix  can  became  quite 
different.  In  particular,  for  the  case  /  ->  °°  the  coefficients  in  the  two  first  columns  go  to  zero  while 
the  other  coefficients  go  to  infinity.  In  such  a  situation  numerical  calculations  on  a  computer  can 
become  unstable  before  achieving  the  necessary  accuracy. 

To  avoid  this  situation,  we  redefine  the  unknown  variables  a/v)  ,  b,(v)  as  new  ones  jt,(v) ,  y/v)  by 
the  formulas 


a?'  =  Af'tfia)  xl2> 


*/ll)  =  Mrii  )yl'> 

b\a)  =  V‘(t\2)  yp 


(3.3) 


Now,  if  we  let  the  unknown  column  matrix  be 

X,  =  (3.4) 

the  system  of  equation  (2.20)  will  be 

W,  X,  =  F,  (3.5) 

where  the  column  matrix  Ft  is  ihe  s.une  as  in  equation  (2.22),  and  VV/  can  be  obtained  from  equation 
(2.23)  by  multiplying  its  columns  by  h,(£,)  ,  ^(q,)  .  ftf‘(£ 2)  •  hr'(r\2)  ,  respectively.  Now  the 
coefficients  of  the  matrix  W/  have  values  of  approximately  the  same  order  (they  differ  from  each  other 
no  more  than  by  a  coefficient  of  order  / )  and  the  numerical  difficulty  mentioned  above  vanishes. 

The  change  in  variable  given  by  equation  (3.3)  is  also  helpful  in  solving  another  numerical  prob¬ 
lem,  that  of  determining  how  much  of  the  infinite  series  of  the  solutions  in  equations  (2.13)  and  (2.14) 
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must  be  retained  in  order  to  achieve  a  given  accuracy.  Note  that  the  unknown  coefficients  now  are  of 
the  same  order  as  the  members  of  Ft 

x,w  -  yr  -  1:(\ 2)  (3.6a) 

and  for  /  >  (;2  we  obtain  die  estimate 


Xl 


(v)  _  v  <v> 
0 


1  - 


/o+2 


'0*7 


-N 


(3.6b) 


where 


e 

/0  =  -e +  N  •  (n>d 


(3.7) 


is  taken  to  be  sufticiendy  large  (/,,  >■  1)  .  Thus,  the  coefficients  x,(v)  and  y/v)  go  to  zero  exponentially 
e 

when  /  >  ~y~  +  N  ;uid  N  -»  °°  ,  Using  the  following  relations  obtained  from  equation  (3.3) 


u,(I)  ],±l(ioaj) 

bi{l)  jt±i(<i$\r) 


y/'to/OliV/ttOlf^Oj 


(  r  <R  ) 


a/2)  hi±t(u>a2r) 
bii2)  hl±]((0$2r) 


-  *i  ^±1(^2  ^VM^) 

=  y/(2>*/  ±i(ti  1  )/*/  (Tiz) 


( r>R  ) 


and  the  "limited"  character  of  the  values 


hti&jit&j) 


it  is  clear  that  the  series  in  equations  (2.13)  and  (2.14)  converge  at  a  rate  that  is  equal  to  or  greater  than 
the  series  in  equation  (2.12).  Therefore,  if  we  take  the  first  / 0  members  of  the  series,  with  1()  specified 
by  equation  (3.7),  the  absolute  error  of  the  calculations  will  not  be  greater  than  a  value  of  order  e~N . 
When  any  of  the  arguments  ,  qv  is  no  more  than  ^2,  the  error  in  the  respective  series  decreases  with 
increasing  r  as  (r/R )"  in  medium  v  =  1  (r  <  R)  and  as  ( Rlr )'°  in  medium  v  =  2  (r  >  R)  .  Thus, 
given  a  error  tolerance  for  the  calculations  of  the  fields  Uv(r  ,0),  we  can  estimate  a  number  l0  for  any 
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frequency  to  and  sum  the  series  in  equations  (2.13)  and  (2.14)  for  /  <  /0.  Test  calculations  showed 
that  V  =  15  in  equation  (3.7)  is  sufficient  to  give  an  accuracy  of  1(T*  or  better. 

As  mentioned  at  the  beginning  of  this  section,  the  system  of  equations  (2.20),  or  equivalently 
(3.5),  are  generally  solved  by  numerical  methods.  However,  analytical  solutions  of  these  systems  are 
possible,  although  laborious,  requiring  the  equivalent  of  inverting  a  matrix  of  rank  four.  This  has  been 
done  and  the  results  are  given  in  Appendix  B.  Having  two  different  solutions  of  the  same  basic  equa¬ 
tions  is  quite  useful  in  checking  the  stability  and  accuracy  of  both  soiutioas.  Thus  the  analytical  solu¬ 
tions  of  Appendix  B  were  used  to  check  the  solutions  that  were  obtained  by  numerically  solving  equa¬ 
tion  (3.5).  This  comparison  showed  that  for  /u  given  by  equation  (3.7)  with  N  =  15  the  numerical  solu¬ 
tion  of  equation  (3.5)  is  suable  and  gives  the  same  result  as  the  analytical  solutions.  For  N  >  15  (accu¬ 
racies  better  than  I0-8)  the  calculations  revealed  that  the  analytical  solution  gave  more  accurate  results. 
It  should  be  mentioned  here  that  all  of  the  calculations  discussed  in  this  paper  were  done  with  double 
precision  arithmetic. 

There  is  another  way  of  estimating  the  accuracy  of  the  calculations  without  considering  the  pro¬ 
perties  of  the  coefficients  Hi  and  ht .  This  approach  hikes  advantage  of  the  observation  made  earlier  that 
there  is  no  coupling  between  harmonics  associated  with  different  values  of  /.  Thus  one  considers  the 
expression  for  the  incident  wave  in  equation  (2.12)  evaluated  on  the  surface  of  the  sphere  (r  =  R)  and 
truncates  the  series  at  a  value  of  /  that  gives  the  desired  accuracy  in  representing  the  incident  wave. 
Then  the  at  and  ht  summed  up  to  this  same  value  of  /  can  be  considered  as  either  the  exact  solution  to 
the  truncated  version  of  the  incident  wave  or  as  the  approximate  solution  to  the  an;ilytical  version  of  the 
incident  wave  having  the  desired  accuracy. 

4.  Scattering  cross  sections 

A  useful  method  of  characterizing  scattering  by  an  object  is  to  calculate  the  energy  of  the  scat¬ 
tered  waves  and  compare  it  to  the  energy  of  the  incident  wave.  Various  forms  of  this  ratio  between  the 
scattered  and  incident  energies  are  called  scattering  cross  sections.  The  energy  of  the  scattered  waves 
can  be  obtained  by  calculating  the  energy  (lux  of  scattered  waves  through  a  surface  S  that  completely 
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surrounds  the  object.  Noting  that  the  energy  flux  per  unit  time  through  a  surface  element  els  having  a 
normal  n  is  given  by  (U  t„  [U])  and  that  the  energy  flux  averaged  over  one  period  is 
to  Im(U  •  tTI*[(J])/2,  then  the  total  energy  flux  per  period  through  the  surface  S  is  given  by 

F  =  yImj(U  OU])  els  (4.1) 

where  (*)  means  the  complex  conjugate.  The  same  reasoning  applied  to  the  incident  P  wave  yields  an 
energy  flux  per  unit  area  of  the  wave  front  of  co2a3(Xi+2)ii)/2. 

Now  we  follow  Truell  et  al.  (1969)  and  define  a  normalized  cross  section  as  the  ratio  of  the  flow 
of  total  energy  carried  outward  by  the  scattered  waves  to  the  rate  of  flow  in  the  incident  wave  through  a 
normal  area  equal  to  the  cross-sectional  area  of  the  object  (geometric  shadow  of  the  object).  In  our 
case  we  let  S  be  a  spherical  surface  of  radius  r„  with  r„  >  R  and  the  normalized  cross  section  is  given 
by 


oN 


Im 


d>  ds 

coa2(X2+2n2)7if?2 


where  can  be  calculated  by  using  either  of  the  formulas 
d>  =  (ll2  •  lf*[U,]) 


(4.2a) 


or 


d>  =  -  (U2  ■  t/fU(>l)  -  (U„  •  t/[U2)) 


(4.2b) 


The  second  expression  for  <t>  is  easily  obtained  from  the  fact  that  the  scatterer  generates  no  additional 
energy,  which  leads  to  the  conditions  that  for  both  the  incident  wave  U„  and  the  total  field  U„+U2  there 
is  zero  net  energy  flux  across  the  surface  S  so  long  as  it  is  entirely  contained  in  medium  v=2  . 


Using  equation  (2  4)  for  the  traction  vector,  the  series  in  equations  (2.11)  and  (2.14)  for  U„  and 
U2,  and  orthogonal  properties  of  the  spherical  vectors,  we  can  obtain  two  different  expressions  for  aN 


On 


4  £<2/  +  1) 

/>> 


} 


+  Id  +  1)Y2 


(4.3a) 


or 


O/v 


-  4 


5>/  +  1) 

1 20 


Re  a/2> 

~w~ 


(4.3b) 


This  result  together  with  the  earlier  noted  fact  that  the  coefficients  for  different  values  of  /  are  not  cou¬ 
pled  means  that  we  must  have 


a. 


,(2) 


+  /(/  +  1)Y23 


bp> 


=  -  Re  a/31 


(4.4) 


for  any  value  of  / .  This  relationship  is  useful  in  verifying  numerical  calculations. 

The  function  from  equation  (4.3)  is  closely  connected  with  the  value  of  the  field  U2(r,0)  of 
scattered  waves  that  is  observed  along  the  positive  z  axis  for  large  r  (r  >  2nVp  /oo)  .  Using  asymptotic 
expressions  for  the  spherical  Hankel  functions  in  equation  (2.14)  we  obtain 


U2(r,0) 


R  YP  +  1)  a/21 
120 


1  +  0(~) 
r 


1  +  0(-) 
r 


(4.5) 


where 


au  =  ij-  Y21  +  »> 


1 2(1 


Comparing  equation  (4.5)  with  (4.3b)  we  have 


aN  =  -4 


Im  A„ 


(4.6) 


which  is  iinalogous  to  ;ui  optical  theorem  for  elastodynamics.  This  expression  in  equation  (4.6)  is  a 
useful  relationship  between  the  amplitude  of  the  forward-scattered  field  and  the  total  energy  scattered  in 
all  directioas  by  the  obstacle. 
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5.  Numerical  results 


One  of  the  prim;iry  purposes  of  (his  investigation  was  to  gain  <ui  understanding  of  the  features  of 
waves  scattered  by  a  local  spherical  heterogeneity.  This  problem  can  be  studied  with  approximate 
methods,  such  as  ray  theory,  but  it  is  not  always  clear  if  the  information  provided  about  amplitudes  and 
waveforms  is  accurate.  With  the  results  presented  earlier  in  this  paper,  it  is  possible  to  calculate  these 
scattered  waves  with  the  assurance  that  the  results  are  complete  and  accurate.  Unfortunately,  the  results 
such  as  equations  (2.13)  and  (2.14)  are  not  very  revealing  in  the  form  of  an  analytical  series  and  must 
first  be  converted  to  some  graphical  form,  such  as  synthetic  seismograms  or  scattering  cross  sections, 
before  the  physical  properties  of  the  scattered  waves  become  apparent. 

We  are  interested  in  knowing  what  types  of  waves  are  scattered  by  a  sphere  and  particularly 
interested  in  understanding  where  the  scattered  waves  with  the  largest  amplitudes  can  be  found  in  the 
region  surrounding  the  sphere.  It  is  possible  to  argue  that  the  main  features  of  such  waves  that  are  scat¬ 
tered  by  a  sphere  will  be  approximately  the  same  for  a  wide  class  of  smoothly-shaped  three-dimensional 
heterogeneities,  and  this  broadens  the  utility  of  the  results  presented  in  this  paper  for  providing  insight 
into  realistic  seismological  problems. 

First  consider  the  matter  of  what  types  of  waves  are  scattered  by  the  sphere  and  their  relative 
amplitudes.  These  effects  are  best  observed  in  the  time  domain,  which  requires  that  the  frequency 
domain  solutions  obtained  in  this  paper  be  transformed  to  the  time  domain.  This  was  done  by  applying 
a  numerical  Fourier  transform  to  the  solutions  and  by  assuming  a  broad  hand  pulse  as  the  form  of  the 
input  wave.  This  source  pulse  had  constant  amplitudes  in  the  frequency  band  between  0.00  and  64.0 
Hz  with  a  sampling  interval  of  0  25  Hz,  which  means  that  the  pulse  contained  wavelengths  that  ranged 
from  more  than  one  order  of  magnitude  smaller  than  the  radius  of  the  sphere  to  more  than  an  order  of 
magnitude  larger.  Synthetic  seismograms  were  calculated  for  the  positions  that  are  shown  in  Figure  1, 
which  lie  on  a  line  offset  a  distance  z  from  the  center  of  the  sphere  in  the  direction  of  forward  scatter¬ 
ing  and  extending  from  the  center  of  the  shadow  into  the  fully  illuminated  zone.  Because  of  the  sym¬ 
metry  of  the  problem,  there  are  only  two  non-zero  components  to  the  solutions,  U„  =  Ux(x,0,z,t)  and 
U =  Uz(x,0,z,t). 
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The  solutions  presented  in  this  paper  can  be  calculated  for  arbitrary  elastic  properties  inside  and 
outside  the  sphere,  including  the  case  of  a  fluid  inside  the  sphere.  Here  we  present  results  for  two 
models,  one  representing  a  low-velocity  scatterer  and  the  other  representing  a  high-velocity  scatterer. 
The  elastic  parameters  for  these  two  models  are  as  follows: 


model  1 

- 

Vp0*  =  4.5  km/s  , 

Vs(l>  =  2.6  km/s  , 

P!  =  2.3  gm/cm 

Vp(2)  =  6.0  km  Is  , 

V/2)  =  3.5  km/s  , 

p2  =  2.7  gm/cm 

model  2 

- 

Vpil)  =  7.5  km/s  , 

Vjl)  =  4.4  km/s  , 

p,  =  3.1  gm/cm 

Vp(2)  =  6.0  km/s  , 

V}2)  =  3.5  km/s  , 

p2  =  2.7  gm/cm 

Note  that  there  is  only  one  physical  dimension  in  this  problem,  the  radius  of  the  sphere  R .  Thus  all  of 
the  other  parameters  in  the  problem  can  be  scaled  with  respect  to  this  parameter,  which  is  given  a  value 
of  unity  in  the  results  (hat  follow.  The  velocities  and  frequencies  scale  with  R  in  the  sense  (hat  the 
results  are  invariant  so  long  as  the  expressions  of  equation  (2.19)  and  the  ratio  r/R  remain  constant. 

Synthetic  seismograms  are  shown  for  three  different  values  of  the  offset  distance  from  the  sphere 
z  in  Figures  2,  3,  and  4.  The  results  for  the  low  velocity  scatterer  (model  1)  are  shown  in  the  upj>er 
parts  of  these  figures  and  those  for  a  high  velocity  scatterer  (model  2)  are  shown  in  the  lower  parts.  In 
all  of  these  figures  a  distance  of  1  km  along  the  x  axis  represents  the  edge  of  the  geometrical  shadow, 
since  the  value  R  =  1  km  was  used  for  the  radius  of  the  sphere 

It  is  useful  to  associate  geometrical  ray  paths  with  the  primary  features  of  the  synthetic  seismo- 
grtuns.  Within  the  forwitrd  shadow  (here  are  three  main  types  of  scattered  waves.  These  ;ue: 

1)  A  compressional  refracted  wave  P2P\P7  that  goes  through  the  sphere.  2)  A  shear  refracted  wave 
that  goes  through  the  sphere  as  P2P \S2  for  a  low-velocity  sphere  and  as  P2.V1.S2  for  a  high-velocity 
sphere.  3)  A  compressional  diffracted  wave  P2P7P7  that  goes  around  the  surface  of  the  sphere.  The 
relative  contributions  of  these  different  waves  to  the  total  seismogram  depends  upon  the  offset  distance 
z  from  the  center  of  the  sphere 

For  nc;ir  offset  distances  :  =  2 R  (Figure  2).  the  wave  P:P  tP2  dominates  the  lost  arrivals  within 
the  shadow  on  the  z  component.  For  a  low  velocity  sphere  (model  1)  it  has  a  reversed  polarity  and 
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arrives  alter  the  wave  diffracted  around  the  outside  of  die  sphere.  Note  that  there  is  a  strong  focusing 
along  die  axis  for  die  low-velocity  sphere  hut  not  the  high-velocity  sphere.  The  diffracted  wave 
P 2P  2P 2  has  a  smaller  tunplilude  within  the  deep  shadow  at  this  offset  distance.  This  low-frequency 
wave  loses  amplitude  exponentially  along  its  ray  path  around  the  sphere  and  is  delayed  in  dme  with 
respect  to  the  undisturbed  incident  wave.  The  shear  waves  P2P  \S2  (low-velocity  case)  and  P2S  \S2 
(high-velocity  case)  are  strong  on  the  x  component  and  have  a  compound  wave  structure,  involving  a 
caustic  and  two  diffracted  waves.  The  caustic  is  due  to  the  fact  that  the  sphere  is  a  low-velocity  zone 
for  both  of  these  waves,  and  extending  beyond  this  caustic  are  low-lrequency  diffracted  waves  that 
attenuate  rapidly  with  distance.  For  the  low-velocity  sphere  this  diffraction  from  the  caustic  begins  to 
interfere  with  the  P2P\P2  wave  at  this  offset  distance.  Behind  this  caustic  and  its  diffraction  is  the 
wave  that  has  been  refracted  through  the  opposite  side  of  the  sphere  and  thus  arrives  at  a  later  time. 
Along  this  branch  of  the  wave  there  is  a  transition  from  a  geometrical  ray  arrival  to  a  diffracted  wave 
that  continues  on  to  greater  distances,  although  the  point  of  this  transition  is  not  obvious  on  the  seismo¬ 
grams.  Note  that,  because  of  focusing  effects  and  near  field  terms  that  are  important  at  this  offset  dis¬ 
tance,  these  waves  have  significant  amplitudes  on  the  z  component  of  motion  near  the  center  of  the  sha¬ 
dow.  In  fact,  at  x  =  0,  we  have  the  unusual  situation  of  an  S  wave  that  appears  only  on  the  z  com¬ 
ponent. 

For  medium  offset  distances  ;  =  4R  (Figure  3),  all  three  of  the  waves  P2P \P2,  P2P2P2,  and 
P2S  {S2  have  comparable  amplitudes.  The  waves  P2P\P2  and  P2P2P2  arrive  closely  in  time  and  inter¬ 
fere  with  each  other.  Ihe  P2P  ,.S\  and  P2SlS2  waves  are  now  separated  in  time  from  the  other  phases 
and  their  compound  nature  is  clearly  evident. 

For  the  larger  offset  distances  :  =  8/?  (Figure  4),  the  wave  P2P  2P  2  is  dominant  on  the  z  com¬ 
ponents,  as  the  shadow  is  not  nearly  so  deep  as  at  the  smaller  offset  distances.  The  wave  P2P\P2 
which  has  passed  through  the  sphere  is  relatively  small  compared  to  the  wave  which  has  diffracted 
around  it  On  the  x  components  the  waves  P2P\S2  and  P2S\S2  are  the  dominant  arrivals  and  still 
display  their  compound  wave  forms.  It  is  worth  commenting  on  the  reason  why  the  P2P\S2  wave  is 
the  dominant  shear  wave  arrival  for  the  low-velocity  case,  while  P2StS2  is  dominant  for  the  high- 
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velocity  case.  These  are  the  S  waves  that  see  the  smallest  change  in  material  properties  upon  both 
entering  and  leaving  the  sphere,  and  thus  the  amount  of  energy  lost  to  reflection  and  also  the  amount  of 
strong  focusing  is  minimized.  In  contrast,  the  wave  P  S  iS2  for  the  low-velocity  case  encounters  a  large 
change  in  velocity  upon  entering  the  sphere,  which  leads  to  more  energy  lost  to  reflected  waves  and 
also  causes  focusing  so  strong  that  the  focal  point  lies  within  the  sphere.  Likewise,  the  wave  P2P\S2 
for  the  high-velocity  case  encounters  a  large  change  in  velocity  upon  leaving  the  sphere. 

It  should  be  emphasized  that,  regardless  of  the  distance  from  the  scatterer,  there  is  a  significant 
amount  of  motion  on  the  x  component.  Recall  that  the  displacement  on  this  component  would  he  zero 
if  the  sphere  were  not  present,  so  we  see  that  scattering  is  an  effective  means  of  transferring  motion 
from  one  component  of  motion  to  an  orthogonal  component  in  the  direction  of  forward  propagation.  At 
the  larger  offset  distances  and  for  both  the  low-velocity  case  and  the  high-velocity  case  the  shear  waves 
are  the  dominant  contributor  to  the  x  -  component  seismograms. 

A  slightly  different  presentation  of  some  of  the  data  from  Figures  2,  3,  and  4  is  shown  in  Figure 
5.  Here  the  seismograms  for  the  z  component  at  the  center  of  the  shadow  are  shown  as  a  function  of 
the  offset  distance  z.  They  show  quite  clearly  how  the  dominant  wave  type  in  the  shadow  changes 
from  P 2P  \P 2  at  small  z  to  PiPjPz  at  large  z  for  the  cases  of  both  the  low  velocity  inclusion  and  the 
high  velocity  inclusion.  Note  that  for  the  low  velocity  inclusion  the  P2P2P2  wave  is  the  first  arrival, 
whereas  for  the  high  velocity  sphere  it  is  the  P2P\P2  wave  that  arrives  first.  This  figure  also  illustrates 
how  the  low  velocity  sphere  is  much  more  effective  in  focusing  energy  into  the  shadow  than  is  the  high 
velocity  sphere. 

These  results  show  that,  depending  upon  the  distance  of  the  observation  point  from  the  scattering 
object,  the  dominant  part  of  the  seismic  field  in  the  shadow  zone  may  be  composed  of  waves  having  a 
fundamentally  different  nature.  Moreover,  in  the  case  of  low  frequencies  it  may  be  difficult  to  separate 
these  different  waves  from  each  other  on  the  basis  of  travel  times.  It  would  appear  that  results  of  this 
type  would  be  applicable  to  various  types  of  seismic  tomographic  methods  For  instance,  the  resulLs  in 
Figure  show  quite  clearly  that  at  observation  points  distiuit  from  the  inclusion  more  than  a  few  times 
its  dimension,  the  dominant  wave  has  been  diffracted  around  the  obstacle  rather  thru)  passing  through  it. 
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which  differs  from  the  common  assumption  made  in  seismic  transmission  tomography.  Because  of  this, 
the  travel  time  anomaly  will  decrease  with  distance  from  the  inclusion.  Obviously,  when  solving 
inverse  problems  to  obtain  estimates  of  the  properties  of  the  scattering  object,  it  is  important  to  have  a 
proper  understanding  of  what  types  of  waves  tire  dominating  the  seismograms  in  order  that  the  correct 
algorithms  can  be  applied. 


So  far  we  have  considered  only  forward  scattered  waves  and  the  results  have  been  illustrated  with 
broad  band  time  domain  seismograms.  To  obtain  a  better  understanding  of  how  the  scattered  fields 
depend  upon  frequency  tutd  also  to  expand  the  spatial  coverage  to  include  back  scattered  waves,  it  is 
convenient  to  consider  scattering  diagrams.  The  field  U2  from  equation  (2.14)  can  be  represented  in  the 


u2  =  fiye)l  f  +  fiye)l  e  +  fu,<e)l  f  +  fu,(e)]  §  (5.i) 

l  Jr  L  JO  l  JT  l  J  ® 

Then,  using  equation  (2.14)  and  the  definitions  of  the  spherical  vectors,  we  obtain  the  following  expres¬ 
sions  for  the  scattering  diagrams  for  the  P  and  S  fields 


f  j  X"1  /•>,  -i*<  r+1 » d/»,  <4) 

iye)Jf  =  -  2/'/:’(2/  +  l)e  3  P,(cosd) 


-'7«+ Utility  dPi(cosQ) 


[U^(0)JB  =  -  £«/3)(2(  +  l)e  2 


4  dQ 


(  1  m  /j/(r|) 

U,(0)  =  y b,l2> /(/+l)(2/  +  l)e  2  -  P, (cos 9) 

>  rr! 


[u>(0)]«  =  5>'2,(2/+1) 


dfy(Tl)  M*l)  1  dPt  (cos  6) 


where  we  have  assumed  that  the  parameters 


,  _  to  r 

_  0(2)  ' 
w> 


;ire  evaluated  for  any  r>R 


19 


Scattering  diagrams  for  six  different  ratios  of  R/X  that  range  from  less  than  0.1  to  greater  than  3.0 
and  for  both  models  are  shown  in  Figure  6,  where  the  radius  of  the  observation  point  r  is  taken  large 
enough  that  near-field  terms  are  small  and  the  scattered  P  and  S  waves  have  their  natural  polarizations. 
A  number  of  interesting  results  emerge  from  the  study  of  these  figures.  First,  at  low  frequencies  more 
of  the  incident  P  field  is  converted  to  scattered  S  fields  than  to  scattered  P  fields,  while  at  high  frequen¬ 
cies  just  the  opposite  occurs.  Second,  at  low  frequencies  the  portions  of  the  incident  field  which  are 
forward  scattered  and  back  scattered  are  comparable,  whereas  at  high  frequencies  most  of  the  scattered 
field  lies  in  the  forward  direction.  At  the  highest  frequencies  we  approach  the  case  of  only  generating  a 
forward  scattered  P  field.  1'hird,  in  terms  of  the  shapes  and  amplitudes  of  the  scattering  diagrams,  there 
are  only  minor  differences  between  the  case  of  the  low  velocity  inclusion  and  that  of  the  high  velocity 
inclusion.  This  is  consistent  with  the  results  shown  in  Figures  2-5,  where  the  differences  between  the 
low  and  high  velocity  case  are  most  pronounced  at  small  distances  from  the  inclusion  and  tend  to 
diminish  at  larger  distances.  Obviously,  the  near  field  parts  of  the  solution  are  more  sensitive  to  the 
sign  of  the  velocity  anomaly  of  the  inclusion  th;ui  are  the  far  field  parts. 

Finally,  to  obtain  even  more  detailed  information  about  the  frequency  dependence  of  the  scattered 
fields,  it  is  instructive  to  consider  the  scattering  cross  sections  defined  in  Section  4.  These  scattering 
cross  sections  are  plotted  for  both  the  low  velocity  inclusion  and  the  high  velocity  inclusion  in  Figure  7. 
An  interesting  result  shown  here  is  that  at  low  frequencies  (kpR  <  2)  more  energy  is  scattered  into  the 
S  field  than  into  the  P  field,  with  the  ratio  reaching  a  factor  of  2  or  more  in  some  ranges.  At  high  fre¬ 
quencies  (kr R  greater  than  about  5).  the  scattered  energy  flux  is  primarily  in  the  P  field.  In  this  range 
the  long  large  amplitude  oscillations  are  caused  by  interference  between  P2P]Pi  and  P2PiPz  waves. 
The  short  small  iunplitude  oscillations  that  appear  only  for  the  case  of  the  low  velocity  sphere  are 
caused  by  the  focusing  characteristic  of  a  low-velocity  inhomogeneity,  essentially  the  waves  multiply 
rellcctcd  within  the  sphere  that  ;ire  evident  in  Figure  5 

In  both  the  scattering  diagrams  of  Figure  6  and  the  scattering  cross  sections  of  Figure  7  the 
iunplitudes  of  the  scattered  fields  have  been  displayed,  but  information  about  their  phases  are  not 
shown  It  is  important  to  remember  that  outside  the  sphere  the  total  solution  consists  of  the  sum  of  the 
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incident  field  U„  and  die  scattered  lield  U2,  and  these  two  fields  can  add  either  constructively  or  des¬ 
tructively.  From  a  physical  viewpoint  it  is  clear  that  when  the  primary  field  U„  interacts  with  the  inclu¬ 
sion,  it  loses  that  part  of  its  energy  which  is  converted  to  scattered  waves  and  this  causes  a  change  in 
the  primary  field.  Thus,  the  additional  held  U2  must  include  both  this  change  in  the  primary  field  as 
well  as  the  secondary  scattered  waves.  For  instance,  in  the  deep  shadow  where  llie  total  field  is  approx¬ 
imately  zero,  the  field  U2  must  have  a  value  that  is  close  to  -  U(,  in  order  to  compensate  for  the  pri¬ 
mary  field.  Consequendy,  in  the  scattering  diagrams  of  Figure  6  the  forward  scattered  P  field 
approaches  a  value  of  1  at  high  frequencies,  hut  this  is  mostly  present  in  order  to  cancel  the  primary 
wave  so  that  the  total  field  present  m  this  region  is  actually  quite  small.  In  the  case  of  the  scattering 
cross  sections  of  Figure  7  this  same  phenomenon  causes  the  scattered  P  field  to  approach  an  asymptotic 
value  of  approximately  2  at  high  frequencies.  This  follows  from  the  fact  that  in  this  frequency  range 
the  scattered  field  consists  of  two  basic  terms,  the  energy  necessary  to  cancel  the  primary  field  in  the 
shadow  and  the  energy  which  is  scattered  by  the  sphere,  and  by  conservation  of  energy  these  two  terms 
must  be  equal  to  each  other  and  equal  to  the  incident  energy  which  is  used  to  normalize  the  cross  sec¬ 
tions  It  helps  in  understanding  this  phenomenon  to  consider  the  case  of  a  perfectly  absorbing  sphere, 
in  which  case  there  will  not  be  any  waves  that  are  actually  scattered,  but  the  normalized  scattering  cross 
sections  will  still  have  a  value  of  1. 

6.  Conclusions 

The  results  presented  in  this  paper  can  be  used  in  a  number  of  different  ways.  First,  the  complete 
and  exact  solution  for  the  scattering  of  an  elastic  P  wave  by  a  spherical  inclusion  is  presented  in  a  con¬ 
venient  form  The  various  terms  of  the  solution  can  be  obtained  either  from  the  analytical  solutions  or 
a  stable  numerical  procedure  Furthermore,  a  simple  expression  is  developed  for  estimating  the 
number  of  terms  required  in  the  harmonic  series  solution  The  results  are  appropriate  for  arbitrary  non¬ 
absorbing  spherical  inclusions,  including  fluids,  and  for  all  frequencies. 

The  numerical  results  calculated  for  the  waves  scattered  from  a  spherical  inclusion  display  a  rich 
variety  of  interesting  phenomena,  some  of  which  were  briefly  investigated  in  this  study.  In  addition  to 
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the  wavelength  of  the  incident  wave,  the  results  depend  critically  upon  the  distance  of  the  observation 
point  from  the  scattering  object.  Within  the  shadow  zone  the  waves  which  dominate  the  solution  near 
the  scatlerer  are  quite  different  from  those  that  dominate  at  larger  distances.  The  solutions  for  low 
velocity  inclusions  are  quite  different  from  those  for  high  velocity  inclusions  when  near  the  scatterer, 
but  these  differences  diminish  at  larger  distances.  Also,  the  solutions  demonstrate  that  an  inclusion  can 
be  quite  effective  in  transforming  incident  P  waves  to  scattered  S  waves,  which  may  be  important  in  the 
formation  of  seismic  codas. 

The  results  presented  in  this  paper  should  also  be  applicable,  in  at  least  an  approximate  manner, 
for  problems  involving  scattering  from  heterogeneities  more  complicated  than  a  simple  sphere.  Scatter¬ 
ing  by  a  sphere  serves  as  a  canonical  problem  for  a  general  class  of  objects  with  relatively  simple  and 
smooth  boundaries  It  is  expected  that  many  of  the  scattering  phenomena  associated  with  the  spherical 
inclusion  will  also  apply  to  this  wider  class  of  heterogeneities. 

Finally,  as  mentioned  in  the  introduction,  the  results  of  this  study  can  serve  as  the  starting  point 
for  an  investigation  of  various  approximations  that  are  typically  made  in  the  interpretation  of  scattered 
seismic  waves.  In  assessing  the  validity  of  these  approximations  it  is  important  to  make  comparisons 
with  the  exact  solutions. 

Appendix  A  -  Spherical  vector  system 

The  spherical  vector  system  used  in  this  paper  and  its  application  to  the  problems  of  elastodynam- 
ics  were  developed  by  G.  I  Petrashen.  His  original  papers  were  published  many  years  ago  and  in  Rus¬ 
sian.  and  are  not  readily  available  today  Thus  we  have  included  a  brief  summary  of  the  system  and  its 
main  properties  in  this  appendix  More  information  can  be  tound  in  Petrashen  (1945,  1949)  and  Kor¬ 
neev  and  Petrashen  (19S7). 

In  a  spherical  coordinate  system  (r.9,<J>)  with  unit  vectors  (r,0,$)  the  spherical  vectors  are 
debited  by  the  expressions 

Y/:  =  K,m)  =  r  X  VK,„,<e.0) 
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Y/m  ■  Y,;(9,<t»)  =  (/  +  !)  f  K/m(0,<t>)  -  r  VK/m(0,<)>) 


Y/m  =  Yta(0.O)  =  1  f  nm(0.<t>)  +  '  Vy'/m(0,(J)) 


with  the  usual  spherical  harmonics 


Y,m  (0.0)  =  />/"(mv0)  .  I  >  0  ,  -/<»!</ 


where  P™(x )  are  the  associated  Legendre  functions 


P/"(x)  =  (l  -xO 


A  2 


-  (X2  _  1)' 


(  m  >  0  ) 


d!ii+'n  2/ 

P,m(x)  =  (  -1  )lml  "--r-!?-  P,'m\x)  ,  (  m  <  0  ) 

(/  +  Im  I)! 


In  a  spherical  c<x>rdinate  system  these  vectors  have  the  components 


Min,  -  1  drlm 
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30  sin0  3(f) 
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In  a  Cartesian  coordinate  system  the  components  are 


V  ^  _ 

*/w  -  — 
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The  spherical  vectors  of  the  system  in  equation  (Al)  are  linearly  independent  at  any  point  (0,4>) 


on  a  spherical  surface.  In  the  space  of  vector  functions  f(0,<J>)  defined  on  a  sphere  Q 
0  <  0  <  rc  ,  0  <  $  <  2n  ,  dQ  =  sin0  c/0  dty 

the  vectors  satisfy  the  following  orthogonality  relation 
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Y*r  •y/J’,  ua  = 


Ci, 


(V) 


-2 


Svv/S//,  6„, 


(A4) 


where  the  normalizing  coefficients  are  given  by  the  expressions 
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JTm 

(/  - 

V  4n  /(/  +  1) 

(l  +  Ml)! 

J  1 

(/  -  MI)! 
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For  vector  functions  f(0,<|>)  with  a  finite  norm 
|  Ifl2  dil  m  |  f*  f  dil  <  - 

the  system  of  spherical  vectors  in  equation  (Al)  is  complete  in  the  sense  of  convergence  in  the  mean 
for  a  generalized  Fourier  series  expansion  of  f(0,cj>) 

-  / 

««.<{»  =  X  X  X  Y/:>(e,<t»  (A6> 

v=0,+,-  /=()  m  -~l 


where 
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(V) 


Y*,r’  ■  f  da 


A7) 


If  the  class  of  functions  f(0,<j>)  possesses  continuous  first  derivatives  with  respect  to  0  and  <t>  on  the 
sphere  Q,  then  the  scries  in  equation  (A6)  convergences  uniformly  and  also  possesses  first  derivatives 
with  respect  to  0  and  <J)  that  agree  with  those  o.  jie  original  function. 
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The  main  feature  of  the  spherical  vectors  of  equation  (Al)  is  contained  in  their  complete  utiliza¬ 


tion  of  central  symmetry  when  acted  upon  by  any  differential  operator  that  is  invariant  to  rotation. 
Consider  an  arbitrary  differential  operator  L  (r  ,8,<j>)  which  is  invariant  to  rotation  of  the  spherical  coor¬ 
dinate  system  (r, 0,4>)  .  By  invariant  to  rotation  we  mean  that  the  operator  L  commutes  with  any  other 
operator  R  that  represents  a  rotation  in  three-dimensional  space.  Then,  for  an  arbitrary  function  \|/(r), 
possessing  sufficient  derivatives,  we  have  the  equality 


L(rM)  [v<r)  Yfifte.*)  =  Af£>y(r)  Y,Lv,<0,$)  ,  v  =  (0,+,-) 


(A8) 


where  M[%\r)  is  an  operator  that  acts  on  the  radial  function  y(r)  only.  This  property  together  with  the 
orthogonality  relation  in  equation  (A4)  means  that  for  a  differential  equation  of  the  form 


L(r,0,<t>)  U  =  0 


substituting  an  expansion  of  the  form 

u  =  X  Yfifte,*) 

vJsn 

converts  the  original  system  into  a  set  of  independent  equations 

AC’  ¥/*>)  =  » 

that  only  involve  the  radial  functions  y ^'(r ). 

In  many  practical  applications,  it  is  necessary  to  know  how  the  vector 

Ufa,  =  y L(r)  Y +  y,;,(r)  \*m  +y fm(r)  Y fm  ,  (A9) 


responds  to  the  following  differential  operators  : 
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Finally  we  note  that  in  the  case  of  boundary  conditions  with  cylindrical  symmetry,  it  is  con¬ 
venient  to  use  a  cylindrical  vector  system 

Y°  =  Y  / 

*  m  1  m  ** 

Ym+  =  Ym  p  -  ,  Ym  +  (A  14) 

v  -  =  Ym  p  +  /  v,„  + 

where  Fm  =  e'"'*  and  (£,$,5)  are  unit  vectors  of  a  cylindrical  coordinate  system  )  .  For  the  case 
of  cylindrical  symmetry  the  vector  system  of  equation  (A  14)  possesses  the  main  features  of  the  system 
of  “phericul  vectors  of  equation  (Al). 


Appendix  B 

The  analytical  solutions  for  the  unknown  coeliicients  in  equation  (2.20)  are 
A(vi 
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The  following  relation  exists  between  the  terms  in  equations  (B-3)  -  (B-5). 
Afc2’  =  f  Im  A^u  Aal)*  j 


An  expression  for  A, l2'  can  be  derived  from  equation  (B.2)  by  substituting  for  functions 
h*(^2)  (k  =  /-l. /,/  +  !)  the  corresponding  functions  -y*(^2)  .  The  following  quantities  were  used  in 
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equations  (B.2)  -  (B.5) 
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0 

Figure  1.  The  geometry  of  the  problem.  The  sphere  of  radius  R  has  material  properties  denoted  by 
v  =  1,  while  the  material  properties  of  the  surrounding  medium  are  denoted  by  v  =  2.  A  plane  P  wave 
is  incident  from  along  the  negative  z  axis.  The  transmitted  wave  field  is  observed  as  a  function  of  the 
x  coordinate  along  a  line  that  is  offset  a  distance  z  from  the  center  of  the  sphere 
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The  upper  panel  is  for  model  1  of  a  low  velcxity  inclusion,  and  ihe  lower  panel  is  for  model  2  of  a 
high  velocity  inclusion 
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Figure  6a.  Scattering  diagrams  tor  various  values  of  the  parameter  kpR  =  (aR/Vp .  In  each  panel  the 
left  hand  figure  represents  the  scattered  P  field,  while  the  right  hand  figure  represents  the  scattered  S 
field.  The  top  two  panels  are  for  the  case  kpR  =  0.5  and  the  bottom  two  panels  are  for  the  case 
kpR  =  1.0.  The  two  panels  on  the  left  are  for  model  1  of  a  low  velocity  inclusion,  and  the  two  panels 
on  the  right  are  for  model  2  of  a  high  velocity  inclusion. 
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Figure  7.  Normalized  scattering  cross  seclitms  as  a  (unction  of  the  parameter  kpR  =  a jR/V,,  (he  lop 
two  panels  are  for  model  I  of  a  low  velocity  inclusion,  and  the  bottom  two  panels  are  for  model  2  of  a 
high  velocity  inclusion.  I  he  panels  on  the  right  tire  expanded  versions  of  those  on  the  left  lor  small 
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Scattering  of  elastic  waves  by  a  spherical  inclusion  -  2.  Limitations  of 

asymptotic  solutions 
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Department  of  Geology  and  Geophysics,  University  of  California, 
and  Center  for  Computational  Seismology,  Lawrence  Berkeley  Laboratory, 
Berkeley,  California  94720 

Starting  with  the  exact  solution  for  the  scattering  of  a  plane  P  wave  by  a  homo¬ 
geneous  spherical  inclusion,  various  types  of  approximate  solutions  are  developed  and 
discussed.  The  standard  Rayleigh  and  Mie  approximations  are  extended  to  the  case  of 
inclusions  having  arbitrary  contrasts  in  material  properties.  For  the  low  contrast  case, 
solutions  are  developed  which  are  valid  over  a  wide  frequency  range.  Several  aspects  of 
these  solutions  are  discussed,  including  the  importance  of  near-held  terms  and  the  rela¬ 
tive  strength  of  the  scattered  P  and  S  fields.  The  various  types  of  approximate  solutions 
are  compared  with  each  other  and  with  the  exact  solution  by  calculating  and  displaying 
their  normalized  scattering  cross  sections. 

Key  words:  elastic  waves,  sphere,  scattering,  diffraction 


1.  Introduction 

Scattering  of  seismic  waves  plays  a  significant  role  in  seismic  wave  propagation  in  the  earth.  Such 
phenomena  as  coda  wave  excitation,  attenuation,  precursors,  and  irregular  arrival  times  are  often  assumed 
to  be  caused  by  scattering  from  inhomogeneities  of  the  earth’s  crust  and  upper  mantle.  The  scattering 
mechanism  also  plays  a  role  in  some  of  the  algorithms  for  migration  and  diffraction  tomography  which  are 
used  in  exploration  seismology. 

F.xact  solutions  for  scattering  problems  are  known  for  only  a  few  types  of  obstacles,  and  even  for 
these  cases  the  calculations  are  not  easily  implemented  because  of  their  complicated  mathematical 
representations  The  basic  approach  to  the  scattering  problem  which  has  been  developed  in  such  disciplines 
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as  optics,  acoustics,  and  quantum  mechanics  makes  extensive  use  of  asymptotic  solutions  to  explain  the 
main  physical  features  of  the  scattered  fields.  In  these  disciplines,  where  the  critical  parameters  of  the 
scattering  obstacle  and  the  observation  system  are  usually  well  known  and  controllable,  experiments  can  be 
arranged  so  that  the  conditions  for  the  asymptotic  solutions  are  strongly  satisfied,  and  thus  these  solutions 
give  a  satisfactory  representation  of  the  observational  results.  However,  in  the  case  of  scattered  seismic 
waves  the  situation  is  more  complicated.  Generally  there  is  a  wide  variety  of  shapes,  locations,  boundaries, 
and  material  properties  in  the  local  heterogeneity  of  the  earth's  crust  and  upper  mantle  which  represent  the 
scattering  obstacles  for  the  seismic  problem.  Additional  complications  arise  from  the  fact  that  these  local 
heterogeneities  are  measured  with  respect  to  a  surrounding  material  which  is  itself  typically  inhomogene¬ 
ous,  but  on  a  larger  distance  scale.  In  most  cases  the  parameters  of  the  scattering  obstacles  are  either 
unknown  or  poorly  known.  Furthermore,  parameters  of  the  observation  system,  such  as  the  locations  of  the 
sources  and  receivers  and  the  wavelengths  of  the  incident  waves,  are  typically  difficult  to  control.  There¬ 
fore,  a  common  situation  in  seismic  problems  is  that  it  is  unknown  whether  the  conditions  necessary  for  the 
application  of  asymptotic  solutions  are  actually  satisfied. 

Establishing  proper  conditions  for  the  use  of  different  asymptotic  approaches  in  the  treatment  of  scat¬ 
tered  seismic  waves  is  not  a  simple  problem  because  of  a  number  of  reasons: 

1.  The  broad  range  in  the  type  and  size  of  scattering  obstacles  which  have  to  be  considered  makes  it 
difficult  to  find  an  asymptotic  approach  which  is  applicable  over  the  entire  range. 

2.  The  fact  that  several  different  parameters  are  involved,  some  of  which  are  poorly  known,  makes  it 
difficult  to  ascertain  whether  an  asymptotic  solution  is  appropriate  or  not. 

3.  It  is  often  necessary  to  apply  a  combination  of  asymptotic  solutions  and  some  of  these  may  be 
incompatible.  For  example,  for  the  far  field  (  r  ->  « )  solution  in  the  case  of  the  Rayleigh  approxi¬ 
mation  ( to  — >  0 )  there  is  the  possibility  of  a  contradiction,  because  the  parameter  kr  =  c or  A'  is 
assumed  to  be  large. 

4.  The  absence  of  a  set  of  canonical  problems  for  which  exact  solutions  are  known  makes  it  impossi¬ 
ble  to  check  the  asymptotic  solutions  by  comparing  them  with  the  exact  solutions. 
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In  the  present  paper  we  eonsider  various  asymptotic  representations  of  the  scattered  fields  formed 
when  a  plane  P-wave  is  incident  on  an  elastic  spherical  inclusion.  An  exact  solution  for  this  problem  and 
methods  of  calculating  it  were  developed  in  a  previous  paper  (Korneev  and  Johnson,  1992).  This  exact 
solution  can  be  used  as  the  starting  point  for  some  of  the  asymptotic  solutions.  It  can  also  be  compared 
with  the  asymptotic  solutions  as  a  means  of  investigating  their  validity. 

2.  Solutions  for  small  inclusions 

We  begin  with  a  brief  review  of  some  of  the  approximate  scattering  solutions  that  are  commonly 
used  in  seismology.  Consider  an  elastic  medium  with  an  inclusion  at  the  center  of  joint  Cartesian  (x  ,y  ,z  ) 
and  spherical  (r.Q.tyj  coordinate  systems.  The  parameters  within  the  inclusion,  denoted  by  the  index 
(v  =  1 ),  are  allowed  to  be  functions  of  the  coordinates 

X,  =  X,(jt.v,z)  ,  p,  =  PiU  ,y.z)  .  p,  =  p,(*,y,z)  (1) 

The  volume  of  the  inclusion  will  be  denoted  by  V.  The  surrounding  media,  denoted  by  the  index  v  =  2 
and  described  by  the  constant  parameters 

X2  =  X  =  const.  ,  p2  =  P  =  const.  ,  p2  =  p  =  const.  (2) 

is  assumed  to  be  homogeneous.  In  the  equations  that  follow,  material  constants  without  indexes  refer  to  the 
surrounding  media  (v=2). 

Most  of  the  approximate  solutioas  used  in  seismology  involve  the  Bom  approximation,  which  is 
assumed  to  be  valid  in  the  case  of  weak  single  scattering  (Aki  and  Chouet,  1975;  Sato,  1984;  Wu  and  Aki, 
1985a,  1985b).  The  actual  conditions  under  which  the  Bom  approximation  is  valid  are  rather  complicated 
and  involve  both  the  dimensions  of  the  inclusion  and  the  contrast  in  its  material  properties  (Hudson  and 
Heritage.  1981),  but  it  is  commonly  assumed  in  seismology  that  they  are  equivalent  to  (he  condition  of  a 
small  contrast  in  material  properties.  By  this  we  mean  that  the  perturbations  of  the  elastic  parameters 
within  the  inclusion  are  small  in  comparison  with  (hose  in  the  surrounding  medium 

|5X|  I  ^1-^2 1  ,  |6p|  I  Pi- P2 1  ,  |8p|  I P 1  P2 1  ,  ,,, 

X  X2  p  p2  p  p2 
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Throughout  this  paper  we  will  be  considering  an  incident  plane  P  wave  of  the  form 


0 


o 


-  it’ 

e  ' 


z  =  U„  e' 


(4) 


This  wave  is  traveling  along  the  z  axis  in  the  positive  direction.  The  solution  to  the  scattering  problem, 
which  includes  both  the  incident  and  scattered  fields,  will  be  written  as 


0  =  Ur'"  =  (U0  +  U,  +  U>*“  (5) 

where  Up  and  If,  represent  the  scattered  P  and  S  waves,  respectively. 

Our  next  approximation  is  to  consider  the  solution  at  low  frequencies.  Then  only  the  first  members 
of  the  frequency  power  series  are  retained  and  the  solution  is 


5A  Sp  „  25p 

- — —  +  cos0  -  - — tr- cox-0 
A  +  2p  p  A  +  2p 


B  -  -^2-  sin0  +  ^~sin20 

P  4 


(6a) 

(6b) 


and  where  the  following  definitions  have  been  used 


.  ,  :  V  e"*'r  „  ,  2  V  0)  ,  to 

A  =  k;  - - ,  8  =  k;  - - ,  k  =  —  ,  ks  =  — 
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The  above  result  has  used  only  the  first  term  in  a  power  series  of  frequency  and  is  thus  dependent  upon  the 
assumption  that 

M  =  ~  «  1  (9) 

vi> 

where  R  is  the  average  radius  of  the  inclusion.  This  is  generally  known  as  the  Rayleigh-Bom  approxima¬ 
tion.  The  solution  of  equation  (6)  obviously  does  not  depend  upon  the  shape  of  the  inclusion,  but  only 
upon  its  volume.  Thus,  for  the  case  of  Rayleigh  scattering,  all  parts  of  a  homogeneous  inclusion  contribute 
in  proportion  to  their  volume  and  the  scattered  field  is  a  simple  sum  of  these  contributions. 
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When  the  wavelength  of  the  incident  wave  is  comparable  to  the  size  of  the  inclusion  and  equation  (9) 
is  no  longer  valid,  then  a  different  approximation  must  be  used.  In  Mie  scattering  the  phase  differences  of 
the  incident  wave  for  different  parts  of  inclusion  are  Liken  into  account  (Chernov,  1960).  The  solution  is 


U,  =  4 


r[  Bessie _  |M)cos!0 

J  A.  +  2ft  p  X  +  2p 


V^l 


mi  *  Y«I,ill2B 
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S,  =  v-'z  -  V-‘f  S2  =  v-'z  -  v;lt 


•V i  =  |s,j  =  21^7 'siity  s,  =  js2j  =  yjv;2  +  V~2  -  2V/71Vi'icos0 


(10a) 

(10b) 


This  approximation  is  equivalent  to  assuming  that  the  inclusion  is  composed  of  numerous  small  nonin- 
leracting  parts,  each  of  which  causes  a  scattered  (ield  of  the  form  of  equation  (6). 

For  the  case  where  the  scattering  volume  is  a  homogeneous  sphere  of  radius  R ,  the  expressions  in 
equation  (10)  can  be  integrated  to  give 
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where 


q  i  =  (aSxR  and  q2  =  <aS2R 
and  ji(q)  is  the  spherical  Bessel  function  of  order  / 

The  solutions  of  equations  (6)  and  (10)  have  been  used  in  many  publications  for  estimating  seismic 
attenuation,  explaining  the  generation  of  coda  waves,  obtaining  high-frequency  asymptotics,  and  formulat¬ 
ing  seismic  diffraction  tomography  In  most  of  these  applications,  however,  valid  bounds  on  the  relevant 
ptirameters  are  not  established,  and  possible  errors  due  to  the  use  of  asymptotic  solutions  are  not 
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considered.  As  a  result,  large  errors  in  the  estimated  seismic  parameters  may  occur,  or  invalid  solutions 
may  be  produced.  In  a  later  section  we  will  illustrate  these  possibilities  by  comparing  the  approximate 
solutions  of  equations  (6)  and  { 10)  with  the  exact  solution  for  a  homogeneous  elastic  sphere. 


3.  Exact  solution  for  the  sphere 

A  method  of  calculating  the  exact  scattering  solution  for  a  homogeneous  elastic  sphere  was 
presented  in  detail  in  a  previous  paper  (Korneev  &  Johnson,  1992),  so  the  results  will  only  be  summarized 
here.  The  solution  for  the  medium  outside  of  the  sphere  for  the  case  of  an  incident  plane  P  wave  has  the 
form 


U  =  U„  +  Uj 


(12) 


where 
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Here  hk(x)  tire  spherical  llankcl  functions  of  the  second  kind.  Analytical  expressions  for  the  coefficients 
iii  and  bi  are  given  in  appendix  B  of  Korneev  and  Johnson  (1992),  as  well  as  information  about  C.j  spheri¬ 
cal  vector  system 

YJL*Y£,(0.0)  =  rx  VK,„,(0,<|» 


(0,40  =  (/+!)  f  Y,m  (0,0)  -  r  V Y,m  (0,0) 


(14) 


Y fm  =  Y,*  (0.0)  =  /  f  Y,m  (0,0)  +  rVYtm  (0,0) 


The  definitions  for  the  spherical  harmonic  functions  are  the  usual  ones 
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Yim (0.0>  =  t m ♦  PPicus 0)  .  /  >  0  ,  (-/  <  m  <  l) 


Any  cylindrically  symmetric  scattered  tield  (as  we  have  for  the  case  of  an  incident  plane  P  wave)  can  be 
represented  in  the  form  of  equation  (13),  where  the  coefficients  «;  correspond  to  the  compressional  field 
and  the  coefficients  bi  correspond  to  the  shear  field. 

For  the  purpose  of  comparing  different  solutions,  we  will  use  the  normalized  scattering  cross  section 
aN,  which  is  the  ratio  of  the  (low  of  the  total  energy  carried  outward  by  the  scattered  waves  to  the  rate  of 
(low  in  the  incident  wave  through  a  normal  area  equal  to  the  cross-sectional  area  of  the  object  (geometric 
shadow  of  the  object).  In  our  case  this  is 


f  (U2  t/[U2j)  ds 
J  axxi(X3+2p2)rri?2 


o#’  +  of}* 


42/2/  +  l)< 

12  0 


|  — 


I-  I  b,  I 

I  +  /(/  +  1)Y2  1 7-=- 1 


(15) 


where  tf[l)2J  is  the  stress  vector  ol  the  field  U2  on  the  spherical  surface  of  radius  r„  ,  (r(t>R)  .  The  two 
parts  of  the  scattering  cross  section,  oj/’’  and  a#’,  represent  the  P  and  S  fields,  respectively.  It  can  be 
shown  that  the  normalized  scattered  cross  section  aN  is  simply  related  to  the  forward  scattered  wave  in  the 
far  field  (r  — »  «>) 
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through  the  formula 
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This  is  the  elastodynamic  equivalent  of  an  optical  theorem. 


(16) 
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4.  Low  frequency  solution  for  an  elastic  sphere  of  arbitrary  contrast 

Starting  with  the  exact  solution  described  in  the  previous  section,  it  is  possible  to  develop  approxi¬ 
mate  solutions  which  are  more  general  than  the  standard  approximations  given  in  section  2.  First  consider 
a  low  frequency  approximation  but  with  no  restrictions  upon  the  contrast  in  material  properties.  In  this 
case  we  retain  from  the  solution  in  equation  (13)  only  the  terms  that  are  of  lowest  degree  in  frequency, 
which  is  to3  and  appears  only  in  the  first  three  (/ =0,1,2)  coefficients 
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lhen  the  solution  has  the  form 
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where  (he  following  functions  have  been  introduced. 
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These  W  functions  depend  upon  frequency  w  and  radius  r  of  the  observation  point  through  the  parameters 
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The  expression  in  equation  ( 19)  is  a  complete  low  frequency  solution  in  that  it  contains  both  the  near 
field  and  far  field  parts  All  of  the  distance  dependence  in  the  solution  is  contained  in  the  W  functions  of 
equation  (20).  We  obtain  the  far  field  asymptotic  form  of  the  solution  if  the  functions  in  equation  (20) 
satisfy  the  following  conditions 

|  Wf;,(2^)  |=l.  \wuz„)  j=l.  ! W$r{Z„)  j  =  1  ,  jwj9(Z,)  j=l.  \w’n{Z.)  j  =  1  (21a) 

\wh(Z„)  j  =  0,  I  W^(ZP)  1  =  0.  j  W\r(Zs)  J  =  0 ,  |w’2,(Zs)J  =  0  (21b) 

The  functions  in  equation  (20)  are  graphed  in  Figure  1  for  the  case  of  y  =  Vs/Vp  =  Zp/Zi  =  1/V3.  This 
figure  shows  that  for  kpr  <  1  the  far  field  approximation  is  definitely  not  satisfied.  As  a)  increases  the 
situation  improves  and  for  kpr  >  3  equation  (21a)  is  satisfied  to  within  10%.  However,  the  terms  of  equa¬ 
tion  (21b)  converge  to  their  asymptotic  values  much  more  slowly  and.  while  three  of  these  terms  have  con¬ 
verged  to  within  10%  of  zero  for  kpr  >  10,  the  W§e  term  does  not  decrease  to  this  level  until  kpr  >  20. 
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What  this  means  is  that  at  observation  distances  r  that  are  only  a  few  times  greater  than  the  wavelength  tl»e 
amplitudes  of  the  scattered  waves  may  he  consistent  with  the  f;ir  held  approximation  but  their  polarizations 
will  be  considerably  more  complicated  than  simple  P  and  S  waves. 

The  situation  is  slightly  more  complicated  than  indicated  in  Figure  1.  In  that  figure  the  absolute 
values  of  the  W  functions  are  plotted,  but  in  equation  (19)  it  is  clear  that  these  functions  can  combine  either 
constructively  or  destructively,  depending  upon  the  difference  in  material  properties  and  the  azimuth  of  the 
observation  point.  Furthermore,  at  small  distances  and  low  frequencies  it  can  be  misleading  to  consider  the 
individual  terms  of  the  solution,  as  what  appears  to  be  a  large  individual  term  may  be  canceled  out  by 
another  and  not  show  up  in  the  total  solution.  Thus  in  Figure  2  the  effect  of  the  near  held  terms  in  the  scat¬ 
tered  held  is  illustrated  by  showing  the  complete  solutions  in  the  time  domain.  Here  the  P  wave  scattered 
from  a  homogeneous  spherical  inclusion  is  shown  at  one  azimuth  and  a  distance  that  is  6  times  the  radius  of 
the  inclusion,  which  corresponds  to  kf,  r  =  1 .2.  The  near  field  terms  make  a  significant  contribution  at  this 
distance,  affecting  both  the  amplitude  and  the  polarization  of  the  scattered  wave. 

In  the  true  far  held  where  all  of  the  conditions  of  equation  (21)  are  satisfied,  the  scattered  held  has 
the  form 
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lliis  result  represents  Rayleigh  scattering  from  ;ui  obstacle  having  an  arbitniry  contrast  in  material  proper¬ 
ties.  In  the  case  where  the  contrast  is  small  and  the  conditions  of  equation  (3)  are  satisfied,  then  equation 
(22)  reduces  to  (6)  and  we  have  the  combined  approximations  of  low  frequency,  far  held,  and  low  contrast. 
Thus  in  the  progression  of  solutions  in  equations  (13),  (19),  (22),  and  (6)  we  see  the  effects  of  the  Rayleigh 
and  Bom  approx  imations. 

The  Mic  scattering  solution  for  a  homogeneous  sphere  in  equation  (11)  can  be  easily  generalized  to 
an  sphere  of  arbitrary  contrast  to  obtain 
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In  order  to  comp;tre  die  v;irious  approximations  discussed  above,  the  nonn;ili/ed  cross  sections  of 
equation  ( 15)  are  plotted  as  a  function  of  kp R  in  Figure  3.  Such  results  are  shown  for  two  different  cases,  a 
low  velocity  inclusion  on  the  top  and  a  high  velocity  inclusion  on  the  bottom.  The  contrast  in  the  velocities 
and  densities  is  about  20%.  The  four  solutions  shown  in  this  figure  are  the  exact  of  equation  (13),  the 
Bom-Rayleigh  approximation  of  equation  (6),  the  arbitrary  contrast  Rayleigh  of  equation  (22),  and  the 
arbitrary  contrast  Mie  of  equation  (23).  Note  that  the  far  tield  approximation  has  been  made  in  all  of  these 
solutions.  In  order  to  calculate  scattering  cross  sections  from  the  solution  in  equation  (23)  it  was  necessary 
to  convert  this  solution  to  the  form  of  equation  (13)  using  the  orthogonal  properties  of  the  spherical  vectors 
in  equation  (14),  and  then  substitute  the  obtained  coefficients  at  and  b,  into  equation  (15).  It  should  be 
pointed  out  that  scattering  cross  sections  for  the  case  of  Rayleigh  scattering  from  a  spherical  inclusion  of 
arbitrary  contrast  were  also  obtained  by  Ying  and  Truell  (1956)  using  potential  functions. 

The  results  shown  in  Figure  3  help  define  the  frequency  ranges  over  which  the  various  approxima¬ 
tions  are  valid.  First  note  that  the  Bom-Rayleigh  approximation  is  the  least  accurate  of  those  shown  in  this 
figure  and  that  the  sign  of  the  error  is  different  for  low-velocity  and  high-velocity  inclusions.  For  the  case 
of  the  low-velocity  inclusion  this  approximation  falls  below  the  exact  solution  at  low  frequencies  and  the 
error  reaches  appreciable  values  in  the  range  around  kpR  =0.4,  whereas  in  the  case  of  the  high-velocity 
inclusion  the  approximate  solution  falls  above  the  exact  solution  and  the  error  exceeds  100%  for  kpR  >  0.3. 
In  contrast  to  this,  the  ;irbitrary  contrast  Rayleigh  approximation  behaves  about  the  same  for  the  low- 
velocity  and  high-velocity  inclusions  and  remains  reasonably  accurate  for  kpR  <0.5.  A  result  which  is 
very  clear  in  this  figure  is  that  by  far  the  best  approximation  to  the  exact  solution  is  provided  by  the  arbi¬ 
trary  contrast  Mie  solution  This  approximation  is  reasonably  valid  for  kpR  <  6  for  the  low  velocity  inclu¬ 
sion  and  for  kpR  <  4  for  the  high  velocity  inclusion,  which  is  an  order  of  magnitude  broader  range  than  for 
the  other  approximations  Of  course,  these  ranges  of  validity  will  vary  somewhat  with  the  values  of  the 
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material  properties  which  are  used,  but  our  calculations  indicate  that  the  results  shown  in  Figure  3  display 
the  primary  features  present  in  the  general  case. 


5.  Low  contrast  approximation 

Another  type  of  approximate  solution  which  can  be  derived  from  the  exact  solution  of  section  3  is 
one  which  assumes  a  small  contrast  in  material  properties  between  the  inclusion  and  surrounding  media, 
but  places  fewer  restrictions  upon  the  applicable  frequency  range  than  does  the  Rayleigh  approximation. 
For  the  case  in  equation  (3)  of  small  perturbatioas  in  the  elastic  parameters  of  the  inclusion,  the  coefficients 
ut  and  bt  of  equation  (13)  (the  original  expressions  can  be  found  in  Appendix  B  of  Korneev  and  Johnson, 
1992)  can  be  reduced  to  the  following  expressions 
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These  solutions  are  not  completely  general  in  frequency,  because  the  manipulations  of  the  Bessel  functions 
that  were  used  in  obtaining  them  are  only  strictly  valid  when  we  have 
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This  restriction  places  an  upper  limit  upon  the  valid  frequency  range  for  equation  (24),  with  this  limit 
increasing  as  the  contrast  in  material  properties  becomes  smaller  and  smaller.  Below  we  will  show  that  this 
restriction  is  significant  for  the  at  coefficients  but  not  for  the  bi  coef  ficients. 


Given  that  the  low  contrast  solution  of  equation  (24)  is  most  appropriate  for  low  frequencies  that 
satisfy  the  restriction  of  equation  (25),  we  can  proceed  differently  to  obtain  solutions  that  are  more 
appropriate  for  the  higher  frequencies.  Here  we  follow  the  method  described  by  Van  dcr  Hulst  (1957)  for 
the  scattering  of  light  from  large  low  contrast  spheres.  The  basic  procedure  is  to  substitute  Debye 
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asymptotic  expansions  for  the  Des.sei  functions  into  the  exact  analytical  expressions.  Starting  with  the 
exact  analytical  expressions  for  the  coefficients  of  equation  (13),  the  low  contrast  assumption  is  m;tde  and 
then  it  is  possible  to  obtain  the  following  approximate  expression  for  the  coefficients  of  the  scattered  P 


wave. 
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and  the  n,  (^)  tire  spherical  Neuman  functions.  Because  we  are  looking  for  solutions  valid  at  high  frequen 
cies,  we  introduce  the  Debye  asymptotic  expansions 


ji<$)  = 


cos  (tj  -  j) 


n,(^)  =  - 


sin (£J  -  — ) 


where 


/  =  sinT  -  tcost  ,  and  cost  = 


'A 


These  approximations  are  valid  for  /  +  —  <  t,.  Now  equation  (27)  can  be  reduced  to 


hi 


1  -  e 


2. it.f  7-tii  1 1  1  -  ike2'^' 


1  +  ike^'f' 


(28) 


53 


where 
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If  the  low  contrast  conditions  of  equation  (3)  are  strongly  satisfied,  equation  (28)  can  be  further 
simplified  to  yield 
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This  is  identical  '.o  the  solution  that  is  obtained  in  optics.  Substituting  the  a,  of  equation  (29)  into  the  opti- 
c;il  theorem  of  equation  (16)  results  in  a  simple  formula  lor  the  normalized  scattering  cross  section  (Moro- 
chnik,  1983). 

°n  -  °ii’)  =  2  -  —  situ/  +  At  (1  -  cos a)  (31) 

a  a2 

An  etirlier  derivation  of  equation  (31)  can  be  found  in  Van  der  Hulst  (1957),  where  it  is  shown  that  this 
result  can  be  explained  by  interference  of  the  incident  and  refracted  waves  propagating  in  the  forward 
direction.  The  parameter  a  is  just  the  phase  difference  between  these  two  waves  in  the  far  field. 

Two  different  formulas  have  been  developed  above  for  the  coefficients  of  the  scattered  P  waves  in 
the  low  contrast  case.  At  frequencies  low  enough  so  that  equation  (25)  is  satisfied,  equation  (24a)  should 
be  used.  At  higher  frequencies  equation  (28)  should  be  used,  and  when  the  low  contrast  assumption  is 
strongly  satisfied  this  can  be  replaced  by  equation  (29).  The  high  frequency  solution  is  valid  so  long  as 

£  <  /  +  y .  The  relationship  between  these  different  approximate  solutions  for  the  scattered  P  wave  is  illus¬ 
trated  in  Figure  4  for  two  different  cases,  a  low  velocity  inclusion  and  a  high  velocity  inclusion.  In  both 
cases  the  contrast  in  velocities  and  densities  is  about  5%.  This  figure  shows  the  scattering  cross  section  of 
the  exact  solution  for  the  scattered  P  field  of  equation  (13)  and  for  two  approximate  solutions,  the  low- 
frequency  P  field  of  equation  (24)  and  the  high-frequency  P  field  of  equation  (28)  Also  shown  is  the 
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scattering  cross  section  which  is  calculated  using  the  analytical  expression  of  equation  (31).  At  low  fre¬ 
quencies  the  approximate  solution  of  equation  (24a)  agrees  best  with  the  exact  solution  for  the  scattered  P 
field,  with  the  error  growing  to  about  30%  by  the  lime  that  the  condition  of  equation  (25)  is  exceeded 
(kpR  -  3).  The  high  frequency  approximate  solution  in  equation  (28)  as  well  as  equation  (31)  both  sys¬ 
tematically  over  estimate  the  exact  solution  in  this  range.  At  higher  frequencies  the  situation  is  reversed. 
The  low  frequency  approximate  solution  of  equation  (24a)  tracks  the  exact  solution  reasonably  well  up  to 
about  kpR  =  5  and  then  is  completely  wrong  at  higher  frequencies.  The  high  frequency  approximate  solu¬ 
tion  of  equation  (28)  and  equation  (31)  agree  quite  well  with  the  exact  solution  for  kpR  >  5,  showing  the 
same  oscillatory  behavior  at  a  slightly  reduced  amplitude  for  the  low-velocity  inclusion  and  at  a  slightly 
increased  amplitude  for  the  high-velocity  inclusion.  In  this  range  these  two  approximate  results  are  so 
close  to  each  other  that  they  appear  as  one  line  on  the  graph.  This  demonstrates  that  the  very  simple 
analytical  expression  of  equation  (31)  provides  an  excellent  approximation  to  the  scattering  cross  section 
for  the  high-frequency  low-contrast  case. 

For  the  scattered  S  waves  in  the  low  contrast  approximation  we  return  to  equation  (24b)  for  the  bt 
coefficients.  In  this  case  the  condition  of  equation  (25)  dtx\s  not  restrict  the  use  of  equation  (24b)  because 
the  higher  order  h,  coefficients  decrease  sufficiently  rapidly  with  increasing  frequency.  At  high  frequen¬ 
cies  the  S  part  of  the  scattering  cross-section.  <3*  \  oscillates  slightly  about  a  constant  level.  Substituting 
equation  (24b)  into  equation  (15)  for  and  numerically  evaluating  the  infinite  sum  of  the  products  of 
Bessel  functions,  we  obtain  the  high-frequency  asymptotic  estimate 
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Ihus  we  see  that  in  the  limit  of  large  frequency,  the  scattering  cross  section  of  the  S  waves  reduces  to  a 
constant  independent  or  frequency. 

An  example  of  the  various  solutions  for  the  normalized  scattering  cross  sections  for  the  S  field  is 
shown  in  Figure  5  It  is  clear  that  the  low  contrast  solution  calculated  from  equation  (24b)  is  a  good 


approximation  to  the  exact  solution  for  both  cases  of  a  low-velocity  inclusion  and  a  high-velocity  inclusion. 
The  agreement  is  extremely  good  for  kpR  <  1,  but  even  at  higher  frequencies  the  approximate  result  tracks 
the  main  features  of  the  exact  result  and  never  differs  from  it  by  more  than  about  25%.  Similar  to  the  scat¬ 
tered  P  field,  this  approximate  solution  under  estimates  the  exact  solution  in  the  case  of  a  low-velocity 
inclusion  and  over  estimates  it  in  the  case  of  a  high-velocity  inclusion.  Also  shown  on  this  plot  is  the  high 
frequency  asymptotic  approximation  of  equation  (32).  While  this  is  obviously  not  an  appropriate  approxi¬ 
mation  at  the  very  low  frequencies,  it  provides  a  simple  but  fairly  accurate  result  at  higher  frequencies.  It 
is  worth  noting  that  the  examples  shown  in  Figures  4  and  5  actually  represent  a  severe  test  for  the  low- 
contrast  approximations.  While  the  contrasts  in  velocities  and  densities  are  only  5%  in  these  examples,  the 
contrasts  in  elastic  constants  are  about  15%. 

Also  shown  in  Figure  5  for  the  purposes  of  comparison  is  the  exact  solution  for  the  normalized 
scattering  cross  section  of  the  P  field.  The  scattered  P  field  is  clearly  much  larger  than  the  scattered  S  field 
for  high  frequencies,  but  for  kpR  <  l  the  situation  can  be  reversed,  with  the  scattered  S  field  considerably 
larger  than  the  scattered  P  field.  Ter  a  certain  degree  this  type  of  behavior  is  a  function  of  the  contrasts  in 
material  properties  of  the  inclusion,  but  there  is  generally  a  range  of  frequencies  where  more  energy  is  scat¬ 
tered  into  the  S  field  than  into  the  P  field.  This  observation  could  have  implications  for  the  generation  of 
seismic  coda. 

6.  Conclusions 

There  are  two  basic  lessons  to  be  learned  from  this  study  First,  considerable  care  must  be  taken  in 
using  various  approximate  solutions  for  the  seismic  scattering  problem  because  these  solutions  are  valid 
only  over  a  limited  range  of  conditions.  Second,  standard  approximations  that  ;tre  typically  made  in 
seismology  and  associated  with  the  names  of  Bom,  Rayleigh,  and  Mie  can  be  considerably  improved  by 
making  a  few  modifications  in  the  formulae  or  by  using  some  of  the  alternative  formulae  developed  in  this 
study  Because  of  the  need  to  make  comparisons  with  exact  solutions,  this  study  has  concentrated  on  the 
case  of  P  waves  scattered  from  a  sphericid  inclusion.  However,  the  results  are  more  genend  than  this  and 
should  be  applicable  to  scattering  Irom  a  more  extended  class  of  objects  with  simple  and  smooth 
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boundaries. 


The  Bom-Rayleigh  approximation  is  commonly  made  in  seismology  for  the  case  of  low  contrast 
scatterers  and  wavelengths  large  compared  to  the  dimensions  of  the  scatterer.  However,  because  of  the 
effect  of  near  field  terms  in  the  scattering  solution,  this  approximation  is  only  valid  at  distances  which  are 
removed  several  wave'engths  from  the  site  of  the  scattering.  Thus  care  must  be  taken  in  using  this  type  of 
an  approximation  to  explain  scattering  in  the  vicinity  of  seismographic  stations. 

Both  the  usual  Rayleigh  approximation  of  equation  (6)  and  the  Mie  approximation  of  equation  (10) 
are  limited  to  low  frequencies  and  low  contrasts  in  material  properties.  The  low  contrast  limitation  can  be 
considerably  expanded  by  using  solutions  developed  as  a  low  frequency  approximation  for  an  inclusion  of 
arbitrary  contrast.  Thus,  by  substituting  equation  (22)  for  the  Rayleigh  formula  of  equation  (6)  and  by  sub¬ 
stituting  equation  (23)  for  the  Mie  formula  of  equation  ( 1 1)  the  range  of  vtilidity  of  these  solutions  can  be 
extended  to  higher  frequencies.  In  particular,  as  shown  in  Figure  3,  the  arbitrary  contrast  Mie  solution  is 
valid  over  a  frequency  range  that  is  an  order  of  magnitude  greater  than  for  the  other  approximate  solutions. 

In  situations  where  the  contrast  in  material  properties  is  actually  small,  approximate  solutions  have 
been  developed  that  are  valid  over  essentially  the  entire  frequency  range.  In  the  case  of  an  incident  P  wave 
two  different  results  must  be  used  for  the  scattered  P  field,  one  for  the  low  frequencies  and  one  for  the  high 
frequencies.  However,  a  single  result  can  be  used  for  the  scattered  S  field  over  the  entire  frequency  range. 
For  situations  where  the  scattering  cross  section  is  needed,  an  analytical  result  based  on  the  optical  theorem 
provides  a  simple  but  good  approximation  for  the  higher  frequencies.  Finally,  these  results  illustrate  some 
of  the  important  differences  between  acoustic  and  elastic  scattering,  as  there  are  frequencies  where  consid¬ 
erably  more  energy  is  scattered  into  S  waves  than  into  P  waves. 
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Figure  2.  Synthetic  seismograms  that  show  the  effects  of  near  field  terms  in  the  scattered  fields.  A  plane  P 
wave  is  incident  upon  a  high-velocity  sphere  for  which  Vp/Vp,=  Vs/V„  =  1.25  and  Pi/p2=  110.  The 
dominant  angular  frequency  of  the  input  pulse  is  to  =  0.2  Vp/R ,  where  R  is  the  radius  of  the  sphere.  The 
seismograms  are  calculated  at  a  radial  distance  r  =  ft/?  and  an  azimuth  0  =  45" .  The  units  of  the  time  t  are 
R/Vp .  The  solid  line  is  the  complete  solution,  whereas  the  dotted  line  includes  only  the  far  field  parts  of 
the  solution. 
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Figure  3.  Normalized  scattering  cross  sections  as  a  function  of  frequency.  The  upper  two  panels  are  for  a 


low  velocity  inclusion  with  Vp/Vpi  =  V,/V3i  =  0.8  and  p,  /p2  =  0.9,  while  the  lower  panels  are  for  a  high 
velocity  inclusion  with  V,,/Vpi  =  V,/V3i  =  1.3  and  p,  /p2  =  1.2.  The  panels  on  the  right  are  expanded  ver¬ 


sions  of  those  on  the  left  for  small  values  of  the  arguments. 
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Figure  4.  Normalized  scattering  cross  sections  as  a  function  of  frequency  for  the  scattered  I*  field  in  the 


low-contrast  case.  The  upper  two  figures  are  for  a  low  velocity  inclusion  with 
V„/V„  =  V,/VSI  =  p,/p2  =  0.95.  while  the  lower  figures  are  for  a  high  velocity  inclusion  with 
Vp/Vpt  =  Vs/V,  =  pi  /p  1  =  1.05.  The  panels  on  the  right  are  expanded  versions  of  those  on  the  left  for 


small  values  of  the  arguments. 
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Figure  5.  Similar  to  Figure  4  for  the  scattered  S  field  Also  shown  for  comparison  is  the  exact  P  field  from 
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